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I. INTRODUCTION 


Orbital prediction has become an essential science needed for several of the DoD 
satellite programs. Exact satellite ephemerides provide for a more accurate means of 
mission analysis. Put more directly, the more accurate the satellite ephemerides can be 
calculated or predicted, the more accurate the mission analysis becomes. The main 
objective being to pinpoint the satellites current or future position. Since the launch of 
Sputnik in the late fifties, a great deal of effort has been placed in trying to model the 
space environment. In particular, the upper atmosphere, or neutral atmosphere, has been 
of great interest in the study of artificial satellite orbit theory. Since 1957, several 
attempts have been made in modeling the thermosphere in order to aid in satellite mission 
analysis. Several atmospheric models, or satellite drag models, are currently used for 
practical applications such as lifetime estimates, reentry prediction, orbit determination and 
tracking, attitude dynamics, and most recently, mass analysis of a particular satellite. 
Atmospheric drag affects all satellites, in all altitude regimes, from low earth orbits to well 
beyond geosynchronous altitudes. For many satellites, the modeling of atmospheric drag 
is the largest error source in describing the forces acting on the satellite. 

Satellite drag models can be divided into two categories, the empirical models and 
the general circulation models. This thesis will compare three of the more widely used 
empirical models. The models used for the comparison are the Jacchia 60, the Jacchia 71, 
and the MSIS 86 earth atmospheric models. Each of these models is formulated in a 
different manner and are unique for the altitude bands that were tested. The ultimate goal 
being to find out which model is more accurate in terms of orbit prediction for a particular 


satellite program. Since the modeling of atmospheric drag ts the largest error in orbital 


prediction, finding the more accurate model for the satellite program in question will 
greatly improve the predicted satellite ephemirides and hence mission analysis. 

The atmospheric models being used in this comparison all differ with respect to the 
input that each model requires. Currently the Air Force Satellite Tracking Center (STC) 
uses the Jacchia 60 atmospheric model which uses not only date and time as inputs but 
also the measurement of the solar flux, F10.7. The J71 model, considered to be an 
improvement of the Jacchia 60 model, adds the measurement of geomagnetic activity, or 
Ap, to its inputs. The MSIS 86 atmospheric model is formulated in a different manner 
than the Jacchia atmospheric models, and uses both F10.7, and Ap as its inputs. 

The atmospheric models calculate the density and constituency of the atmosphere 
based on the current and predicted or average environmental conditions. The accuracy of 
these models has been calculated to be 80 to 85 percent accurate. This percentage drops 
off substantially as the altitude increases. Accuracy also seems to decrease during periods 
of high solar and geomagnetic activity. At the moment, the sun ts on the downside of the 
11 year solar cycle. This is an advantage for the comparison of the atmospheric models, 
because there were fewer and weaker solar flares and geomagnetic storms to perturb the 


upper atmosphere. 


If. BACKGROUND 


A. ATMOSPHERE 

The earth's atmosphere ts classically divided into four different regions based on 
temperature and pressure gradients. These four regions are the troposphere, stratosphere, 
mesosphere and the thermosphere. The corresponding boundary layers, or upper limits of 
each of the regions are the tropopause, stratopause, mesopause, and thermopause 
respectively. Figure 1 (Akasofu, 1972, p.109), illustrates the breakdown of the earth's 
atmospheric regions. Beyond the thermopause is the region delineated as the exosphere. 
The exosphere is a region of extremely low density and temperature, and is the transition 
region into space. 

Another way in which the earth's atmosphere ts divided, is by the classification of 
two regions known as the homosphere and the heterosphere. The transition boundary 
between the regions ts labeled the turbopause. Once again the region outside the 
heterosphere is labeled the exosphere. Figure 2 (Akasofu, 1972, p.111) illustrates the 
various atmospheric regions as derived by the two defining systems. This figure provides 
a breakdown of altitude versus temperature for the various altitude regimes. It should be 
noted at this time, that the atmospheric region of interest during this study was an altitude 
band between 600 to 800 km. This altitude band 1s encompassed within either the 
thermosphere or exosphere, depending on which classification scheme is being used. For 
the sake of continuity, the altitude band of interest will be considered to be within the 


thermosphere. 


1. Troposphere 
The troposphere is that portion of the atmosphere that extends from the surface 
to roughly 10 to 15 km above the surface. It 1s in equilibrium with the sun-warmed 
surface and is characterized by intense convection and cloud formations. In this region, 
both temperature and density decrease with increasing altitude, with an occasional 


inversion layer. (U.S. Air Force, 1960, p.1-3) 


2. Stratosphere 

The stratosphere is located above the tropopause and extends up to 50 km. This 
region of the atmosphere is extremely important in that it contains the ozone that is 
responsible for the absorption of the extreme ultraviolet (EUV) radiation produced by the 
sun. Due to the absorption of this EUV radiation, the stratosphere has a positive 
temperature gradient. The density, however, still decreases with altitude. One other 
consequence of the absorption of the EUV radiation by the ozone layer is that the EUV 
radiation cannot be measured from the earth's surface. This presents a problem which will 


be addressed later in this paper. 


3. Mesosphere 
The Mesosphere 1s that region of the atmosphere located above the stratopause 
and extends up to 80km. Once again, both the temperature and density are decreasing 
with altitude. The mesosphere is in radiative equilibrium between the ultraviolet ozone 
heating by the upper fringe of the stratosphere, and the infrared ozone and carbon dioxide 


cooling by radiation to space. (U.S. Air Force, 1960, p.1-3) 


4. Thermosphere 
The thermosphere extends from the mesopause to higher altitudes with no 


altitude limit. The thermosphere is characterized by a very rapid increase in temperature 


4. Thermosphere 
The thermosphere extends from the mesopause to higher altitudes with no 

altitude limit. The thermosphere is characterized by a very rapid increase in temperature 
with altitude due to the absorption of the sun's EUV radiation. The temperature increase 
reaches a limiting value known as the exospheric temperature, the average values being 
between ~600 to 1200 K over a solar cycle. (Larson, 1992, p.208) The thermosphere may 
also be heated by geomagnetic activity, which transfers energy from the magnetosphere 
and ionosphere to the thermosphere. The heating of the thermosphere causes an increase 
in the atmospheric density due to the expansion of the atmosphere. Figure 3 (Hess, 1965, 
p. 679), illustrates the variation of temperature versus altitude for the various atmospheric 


regimes. 


5. Homosphere 
The homosphere extends from the surface to approximately 100km. It is 

characterized by its uniform composition and relatively constant molecular weight. The 
composition of this region can be broken down into the following: 78% N, 21% O2, 1% 
Ar and trace amounts of other gases. The uniformity of the region is created due to the 
turbulent mixing of the gas constituents. (Adler, 1993, p.10) The composition, hence the 
uniformity, of the homosphere changes at ~100km altitude due mainly to the dissociation 
of the oxygen molecules. Because the density at this altitude is low, recombination of the 
monatomic oxygen is very infrequent; even more so as altitude increases. The dissociation 


of oxygen causes the molecular weight to decrease substantially. 


6. Heterosphere 
The heterosphere exists from ~ 100km outward, with no altitude limit. The 


region 1s characterized by diffusive equilibrium and significantly varying composition. The 


molecular weight of the atmosphere decreases rapidly, from ~ 29 at 90km to ~ 16 at 
500km. Above the region of oxygen dissociation, nitrogen begins to dissociate. Diffusive 
equilibrium begins to take place, and the lighter molecules and atoms rise to the top of the 
atmosphere. The distribution functions, or scale heights. of each of the constituents of the 
heterosphere are found by equating the pressure gradient of the atmosphere with the 


gravitational force, as described by the ideal gas law, 


p=(2)er (1) 


m 
where p = gas pressure, k = Boltzman constant, m = molecular mass, T = temperature, 


and 0 = density. For a small cross sectional area of thickness dh 


dP = —pgdh (2) 
therefore 
aP = ap, at (3) 
jo) we 
ie T{ ae.) (4) 
mg\ p T 


_ , . kT ae 
By assuming isothermal variation and defining the scale height to be 7 = —, a simple 
mg 


differential equation is obtained described by the following; 


ap __1 4, (5) 
p H 
with the solution being; 
| 
-—h 
p=pe" (©) 


Each molecule will have a different scale height depending upon its mass. This gives rise 


to diffusive equilibrium, in that the density of the varying constituents will decrease at 


certain levels. The diffusion process takes place amongst Ar, O07, N>, O, He, and H 
respectively as altitude increases. (Adler, 1992, p.11) Figure 4 ( Akasofu, 1972, p. 110) 


illustrates the densities of the various atmospheric constituents versus altitude. 


7. Time Dependent Variations 
Several time dependent density variations are present in the earth's atmosphere. 
The density of the earth's atmosphere varies according to the time of day, day of the year, 


and which year it happens to be in the 11 year solar cycle. 


a. Diurnal Variation 

The atmospheric density variation which is dependent upon the time of day is 
called the diurnal variation. The maximum density of the upper atmosphere can be found 
at approximately 1400 local time, with the minimum around 0200 local time. The density 
variation becomes more pronounced with an increase in altitude as can be seen in Figure 5. 
(Hess, 1968, p. 99) The diurnal variation is caused by the alternate heating during the day 
and cooling during the night of the upper atmosphere. Often called the diurnal bulge, the 
density variation occurs mainly over the equator, with an elongation in the north-south 
direction due to the tilt of the ecliptic. The peak of the phenomena occurs at the latitude 
of the sub-solar point. (Jacchia, 1963, p. 983) The temperature change in the upper 
atmosphere during this diurnal correction parallels that of the density, except that the 
temperature change lags behind the density change by approximately two hours. This lag 


seems opposite of what would be expected, and 1s not completely understood. 


b. 27-day Variations 
It has been established that there is a density variation that is caused by the 
27-day solar cycle. This was shown by Jacchia, who determined the correlation between 


actual satellite drag measurements and the solar decimetric flux, or 10.7 cm flux. It was 


found that the density of the atmosphere not only had a daily variation, but also a monthly, 
or 27 day cycle, that corresponded to the 27 day rotation cycle of the sun. (Hess, 1968, 


p.101) 


c. Semi-Annual Variations 

One of the least understood upper atmosphere density variations is the semi- 
annual variation. This density variation is characterized by a pronounced minimum during — 
the June-July time frame, with a secondary minimum occurring in January. The 
maximums occur during the September-October time frame, with a lesser secondary 
maximum occurring during March-Apmil. Several theories exist as to what causes the 
semi-annual variation. The most controversial, is that the variation is an effect of the solar 
wind. Another theory is that the variation is caused by the convective flows from the 


summer pole to the winter pole. The flows would descend at the winter pole, transporting — 





heat to the cooler mesosphere from the higher altitudes. (Hess, 1968, p. 103) 


d. Long Term Variations 

It was found that there was a long-term density variation associated with the 
11 year solar cycle. Measurements of the solar flux, or F10.7, taken over several years 
provide evidence of this effect. As can be seen in Figure 6 (Larson, 1992, p. 209), the sun 
has an 11 year cycle with maximum and minimum values of F10.7. The F10.7 
measurement, explained in detail later, is a measurement of the EUV radiation. During 
periods of high solar activity, the solar flux is high, thus causing an increase in atmospheric 
density due to EUV heating and the resultant expansion. As can be assumed from the 
figure, the density of the upper atmosphere has a corresponding eleven year minimum and 


maximum. 


B. THERMOSPHERIC PROCESSES 

The upper atmosphere or thermosphere undergoes a change in composition and 
density due to several external inputs. The majority of the change is caused in response to 
the activity of the sun. The radiation from the sun, both thermal and ultraviolet, cause 
varying rates of change in the composition and hence the density of the atmosphere. The 
sun also causes changes in density due to the solar wind. The impingement of the solar 
material upon the earth's magnetosphere and upper atmosphere causes several changes in 
the make-up of the thermosphere and subsequently the overall density at altitude. The 
other major contributor to density variation of the earth's atmosphere is geomagnetic 
activity. Geomagnetic storms, although short-lived, cause a significant change in the 


atmospheric composition and density. 


1. Solar EUV Radiation Effects 
The first source of density change that will be addressed is the Extreme 

Ultraviolet Radiation produced by the sun. The sun's EUV radiation is the main cause of 
density variation in the earth's thermosphere. The solar EUV radiation is deposited mainly 
at low latitudes in the summer hemisphere. (Marcos, 1993, p. 2) The circulation and 
structure of the thermosphere at the low and middle latitudes are controlled by the heating 
caused by the absorption of the EUV radiation by the ozone layer in the lower 
thermosphere. The longer wavelength UV and visible radiation reach the lower 
atmosphere and hence heat the earth's surface. (Burns, 1991, p. 3) Most of the EUV 
radiation reaching the earth's atmosphere 1s absorbed at the 300 km level and the energy 
enters the atmosphere through photoionization. The energy that has been absorbed by the 
electrons and ions is passed to the neutral atmosphere by collisional processes. (Burns, 


1991, p.3) 


The solar EUV radiation also imparts a momentum to the neutral gas by the 
creation of pressure gradient forces that drive the neutral winds from the day to night 
regions and from the summer to winter hemisphere. (Burns, 1991, p. 3) Variations in the 
strength of the EUV radiation interacting with the thermosphere lead to changes in the 
composition and density of the neutral atmosphere. During the period of solar minimum, 
the EUV radiation produced by the sun is much less than that being produced during solar } 
maximum. Therefore, the thermospheric temperatures and neutral densities will be lower 
during a solar minimum than during solar max. This fact is illustrated in Figure 7. (Larson, 
1992, p. 209) It can be seen from this figure that the variation in density becomes much 
greater at higher altitudes during periods of high solar activity versus low solar activity. 

Due to the fact that the solar EUV radiation is absorbed by the thermosphere, it 
makes it difficult to obtain a measurement of the EUV impinging on the atmosphere. This 
problem was solved by Jacchia in the early sixties. Jacchia discovered that the intensity 
found at 10.7cm closely corresponded to the amount of solar activity being witnessed. 
Currently the accepted measurement of solar flux is the F10.7 index. As can be seen in 
Figures 6 and 7, the typical value for F10.7 during solar minimum is 75, and has reached 
as high as 290 during the solar maximum period of 1958. A change in solar activity of this 
proportion would mean a change in density by a couple orders of magnitude at the altitude 
regime of interest. 

One of the drawbacks of using F10.7 as a gauge of EUV radiation, is the fact 
that it lies at the other end of the spectrum from the EUV radiation, and is not a direct 
measure of the amount of EUV radiation reaching the thermosphere. Figure 8 (Hess, 
1968, p. 668) illustrates the solar spectrum. It can be seen that the EUV radiation Is at a 
frequency on the left hand side of the peak spectrum, whereas the solar flux measurement, 


F10.7, is on the nght. Because the F10.7 index is not a direct measurement of the amount 





of solar EUV radiation entering the atmosphere, accurate density calculations based on 
solar activity have some inherent error. Several programs have been initiated in order to 
remedy this problem, but currently the F10.7 index is the best indicator of solar flux 


available, and consequently is the index most widely used in atmospheric modeling. 


2. Solar Wind 

Another of the driving forces causing composition and density variation in the 
thermosphere is the solar wind. The solar wind consists of protons and electrons flowing 
outward from the sun's corona. Higher density plasma streams are also ejected from the 
sun during periods of flare and sunspot activity. (Fleagle, 1963, p. 236) The solar wind 
blows the interplanetary magnetic field lines across the polar cap in a direction away from 
the sun. (Burns, 1991, p. 4) This in turn causes a potential drop across the earth's 
magnetic polar cap as the interplanetary magnetic field encounters the earth's magnetic 
field. Field-aligned currents flow down to the ionosphere, closing the circuit, and 
producing an ion-convection pattern at high latitudes. The ions in this convection pattern 
collide with the neutral particles, driving them in a similar convection pattern. (Burns, 
1991, p. 4) These collisions produce heat, which in turn produces a rising motion around 
the auroral zone. The up welling and the convection-driven neutral winds produce a 
composition and density change which spreads from the high latitudes to the lower 
latitudes. 

The convection driven neutral winds also produce another significant heat 
source. Joule heating results from the frictional heating of the plasma as it is dragged 
through the neutral upper atmosphere by the auroral electric field forces. (Marcos, 1993, 
p. 3) This Joule heating is a substantial heat source, but becomes even more prevalent 
during periods of solar flare activity. Flares on the sun cause the solar wind to accelerate, 


driving the interplanetary magnetic field faster across the earth's magnetic field. This 


1] 


Increases the cross-cap potential and the rate of particle precipitation, and can also 


produce a magnetic substorm. (Burns, 1991, p. 5) 


3. Geomagnetic Storms 

Another of the major contributors of composition and density variation are 
geomagnetic storms. Geomagnetic storms are produced by the interaction of the 
interplanetary magnetic field with the earth's magnetic field. When solar flare activity is 
occurring, shock waves in the interplanetary magnetic field are driven into the earth's 
magnetic field causing rapid transients in the earth's magnetic field. This effect is 
monitored by means of the Ap index, or geomagnetic activity index, during periods of 
solar flare activity. The geomagnetic storm presents itself as a world-wide transient 
variation in the earth's magnetic field. The onset, or sudden commencement, of a 
magnetic storm 1s characterized by a rapid increase in the Ap index. Approximately 20 
minutes later, the temperature and density of the auroral zones begins to increase. 

One of the manifestations of geomagnetic storms is the generation of waves 
that propagate from the auroral regions to the lower latitudes. These waves take 
approximately eight hours to reach the lower latitudes. (Alder, 1992, p. 18) A typical 
magnetic storm, illustrated in Figure 9 (Akasofu, 1972, p. 557), lasts from 24 to 48 hours 
or longer. The effects of the geomagnetic storm on the density, last even longer and are 
quite pronounced at the higher altitudes. This is illustrated in Figure 10. (Ratcliffe, 1972, 
p. 35) It can be seen from Figure 10, that during periods of geomagnetic storms, the 
density at the altitude regime of interest increases several fold. This fact coupled with 
intense bombardment of EUV radiation, increases the density by several orders of 


magnitude. 





4. Gravity Waves, Planetary Waves, and Tides 

The final source of composition and density variation to be discussed are 
propagating tides and gravity waves. Atmospheric tides are caused primarily by the 
absorption of ultra-violet radiation by the ozone in the stratosphere, while gravity waves 
are caused by a variety of mechanisms which occur in the troposphere. A couple of these 
mechanisms are the shears associated with cold fronts and winds blowing over mountains. 
(Burns, 1991, p. 3) At lower altitudes, the semi-diurnal tide is the major contributor to 
density variation. This effect becomes less apparent at higher altitudes due to the fact that 
the semi-diurnal tide is dissipated by viscosity and ion drag. At higher altitudes, the 
diurnal tide becomes the driver of density variation and is caused by the absorption of the 
EUV radiation in the thermosphere. Overall, these waves and tides propagate up from the 
lower altitudes affecting the composition and density of the upper altitudes. 

As a consequence of the above mentioned thermospheric processes, it is known 
that the temperature and hence the density of the upper atmosphere vary with the 
following parameters: 


e solar flux (solar cycle and daily component) 


geomagnetic activity 
e local time 

e day of the year 

e latitude 

2 longitude 


e wave structures 


13 


C. PROGRAM DEVELOPMENT 

The propagation program created for the atmospheric model comparison is a 
FORTRAN based program which uses a Runge-Kutta integrator. The propagator uses a 
form of Cowell's Method of orbital prediction. The program uses the position and 
velocity vectors in Cartesian coordinates as its inputs, and integrates over a designated 
time frame to produce the position and velocity vectors, also in Cartesian coordinates. 
Cowell's Method of orbital prediction has been found to be inaccurate due to the build-up 
of round off error. A more accurate method would have been to convert the Cartesian 
coordinates into normalized spherical coordinates in order to minimize the integration 
error build-up. Further research into this conversion is being conducted in order to obtain 
greater accuracy. For the purpose of the comparison, however, the accuracy of the result 
was not in question, but the accuracy of the atmospheric model, and its relative ease of 
use. 

The propagation program consists of a series of subroutines which calculate the 
perturbing forces acting on the satellite. These forces, in the form of accelerations, are 
applied to the satellite's motion resulting in a position and velocity vector reflecting the 
result of the perturbing forces. Figure 11 and Table 1(Milani, 1987, p. 14-15) illustrate 
the various perturbing accelerations and their relative magnitude as a function of altitude. 
In order to obtain an accurate reflection of orbital motion at the altitude band in question, 
it was decided that in addition to the drag force, the earth's geopotential and third body 


attractions would also be included as perturbing forces. 





1. Earth's Geopotential 
As can be seen from Figure 11, the main perturbing force encountered by low 


earth orbiting satellites is the earth's gravity field. If the earth was a perfectly round and 
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smooth planet, then the forces of gravity could easily be modeled by the inverse square 


gravity law shown in equation 7. 


g=-r (7) 
r 


However, since the earth is not perfectly spherical, and not smooth, the gravity field must 
be derived by obtaining the solution to Laplace's equation: 

vV =0 (8) 
where 


y=ec{2 (9) 


volume 


6 being the distance from the satellite to an incremental mass, am, inside the earth, and G 
is the factor of proportionality in Newton's Law of Gravitation. Laplace then derived the 


basic equation for the geopotential shown below: 


Gre p\ 
V=— je: —1| di 10 
; Jd (cos0)(2) m (10) 


By converting to spherical coordinates and applying Rodriguez’ formula, the earth's 


geopotential takes the form of equation 11. 


al +yr(2 ) Pa (sin ¢')(C,,, coum sina) (11) 
n=2 m=0 

In this equation, V 1s the gravitational potential function, GM is the earth's gravitational 

constant, r is the radius vector from the earth's center of mass, a is the semi-major axis of 

the model ellipsoid, n and m designate the degree and order of the coefficients, @ is the 

geocentric latitude, A is the geocentric longitude, C_,, and S_ are the harmonic 


coefficients, and P are the associated Legendre functions. (Ross, 1993, p.5-3) Several 


Earth gravitational models are currently in use in both the civilian and DoD communities. 


15 


The DoD currently uses the WGS-84 geopotential model as its gravitational model. The 
WGS-84 Earth gravitational model (EGM) contains a 180 by 180 matrix of zonal and 
tesseral harmonic coefficients. The WGS-84 EGM is currently being used for several 
DoD satellite programs including the GPS satellite constellation, and is considered to be 
an extremely accurate representation of the earth's geopotential. For most cases, the full 
matrix 1s not needed, and the matrix is truncated down to a 41 by 41 matrix. The 
truncated matrix provides an ample representation. 


Currently, the propagation program contains the first 8 sets of zonal and tesseral 





coefficients. Future iterations of the propagation program will contain the complete 41 by 
j 


41 matrix. At the moment, attempts at configuring the propagation program in the WGS- 
84 coordinate system have failed. The first eight sets of coefficients vary only slightly 
from the first eight zonal terms of the more basic geopotential models and can be used 


without incurring significant errors during propagation. 


2. Drag 
For satellites in low earth orbit, drag is the other major perturbing force that 
must be modeled. Drag affects all satellites in all altitude regimes, but the affect is 
considered to be insignificant at altitudes greater than 1000 km. The following equation 


represents the atmospheric drag acceleration which is placed on an orbiting body. 

dy = Cy pV (12) 
In this equation, a, represents the drag acceleration imparted upon an orbiting satellite, 
C, is the coefficient of drag, A is the cross sectional area of the satellite perpendicular to 
the direction of motion, 77 1s the satellite mass, V is the satellite velocity and p is the local 


atmospheric density. The velocity used in this equation is computed by combining the. 


geocentric velocity of the satellite with the contribution due to earth rotation and the wind 


velocity at altitude. It is important to note that the wind component of the velocity cannot 
be ignored, and can be quite significant at altitudes in excess of 600 km. The neutral wind 
is a Consequence of the activity caused by a geomagnetic disturbance. The change in drag 
can be 5% for every 200m/s of wind velocity. During a period of intense geomagnetic 
activity, the neutral winds can reach velocities in excess of 1 km/s which is equivalent to a 
25% change in density. (Marcos, 1993, p. 6) In order to model this behavior, several 
atmospheric wind models have been developed. Models of the neutral winds begin with 
the efforts of Sissenwine in the late sixties, to the Horizontal Wind Model 90 (HWM 90). 
Currently these wind models are not incorporated into the empirical atmospheric drag 
models, but are incorporated in the general circulation models. 

The coefficient of drag, C,, 1s a difficult quantity to obtain. In order to 
determine the coefficient of drag, it must be determined whether the satellite is in a 
continuum flow, or a free molecular flow. This 1s done by determining the Knudsen 
number. The Knudsen number is the ratio between a typical dimension of the satellite and 
the average mean-free-path of the molecules found in the local atmosphere. When the 
Knudsen number is much less than one, the satellite 1s considered to be in a continuum 
flow. The satellite is considered to be in a free molecular flow when the Knudsen number 
is greater than one. Since the atmospheric density is so low at orbital altitudes, satellites in 
the upper atmosphere are characterized by free-molecular-flow aerodynamics. (Ross, 
1994, p.16) 

The collision of the atmospheric particles with the spacecraft produce the 
atmospheric drag force. The collisions can be classified under three categories: (1) elastic 
or specular reflection, (2) diffuse reflection and (3) absorption and subsequent emission. 
In elastic or specular reflection, the molecule collides with the satellite and is reflected 


away without transferring energy to the satellite. In diffuse reflection, the molecule 


by 


collides with the satellite and transfers a portion of its energy. Energy is also transferred 
to the satellite when molecules are absorbed on impact and later emitted. Since the 
satellite is considered to be in a free-molecular-flow, the coefficient of drag is obtained 
from a statistical mechanical calculation. (Ross, 1994, p. 16) Depending on the size and 


makeup of the satellite, the coefficient of drag can vary from 2.0 to 6.0. Table 2 (Larson, 





1992, p. 207) illustrates the variation of the drag coefficient for various orbiting platforms. 
If the coefficient of drag cannot be determined, a value of 2.2 is used. For the purpose of : 
this thesis, the coefficient of drag will remain constant at a value of 2.2. 

In equation 12, the coefficient of drag, the cross sectional area, and the mass of : 
the satellite all make up what is called the B-factor. 


_ Od 
m 


B (13) 


Inverting equation 13 will result in what is more widely known as the Ballistic coefficient. 
Current practice is to guess what the correct B-factor is for that given day and adjust the | 
B-factor until the correct position and velocity is achieved. This does not seem to be 
orbital prediction, but rather orbital correction. If the size and shape of the satellite are 
known, as well as the mass, then the B-factor should only vary with varying coefficients 
of drag. However, most of the current orbital prediction schemes use the B-factor as a 
catch-all for any other errors in the modeling program. Hence the value of the B-factor is 
nowhere near the actual value. In this comparison, it was decided to keep the B-factor at 
a constant value in order to obtain a more accurate comparison between atmospheric 
models. It must be noted at this point that keeping the B-factor at its actual value is a 
major change, vice using the B-factor as an error catch-all. 

The density for the drag acceleration calculations is of course derived from the 


atmospheric models. In the calculation of the drag acceleration, the density 1s the most 
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difficult to obtain, and frequently the one parameter with the greatest error. Asa rule, the 
most current atmospheric models have a 15 to 20% inaccuracy rate, with this inaccuracy 
increasing as altitude increases. Due to this inaccuracy, the propagation scheme is only as 


accurate as the atmospheric model being used. 


3. Third Body Attractions 

The last of the perturbing forces currently included in the propagation program is 
that of third body attractions. In the case of the satellite program in question, the 
perturbing bodies are the sun and moon. As can be seen in Figure 11, both the sun and the 
moon contribute some small portion of disturbance force to a medium altitude satellite. In 
practice, this disturbance force is usually ignored for the low earth orbiting platforms, but 
in the altitude band in question (600 - 800km), these disturbance forces must not be 
ignored if a truly accurate representation of orbital motion is to be modeled. The equation 
used to find the perturbing acceleration due to the moon is described by equation 14 


below. 





: r r 
v --Sr-n,( 3 -“s| (14) 


r a AE; 
In this equation, r is the radius from the earth to the satellite, r,, is the radius from the 
moon to the satellite, -,, 1s the radius from the earth to the moon, and LL, is the 
gravitational parameter of the moon. (Bate, 1971, p. 389) 
Currently the propagation scheme does not contain several perturbations that are 
Important for accuracy purposes. At the moment, the sun's radiation pressure is ignored 
as well as the earth's albedo effect. Also, relativistic effects are ignored, which must 


eventually be incorporated in order to improve accuracy. 
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D. ATMOSPHERIC MODEL DEVELOPMENT 





As mentioned previously, atmospheric models can be divided into two categories, 
the empirical models, and the general circulation models. The general circulation models | 
are dynamic representations of the earth's atmosphere, and require extensive 
computational time in order to model the atmosphere. The general circulation models that 
have been recently created require Cray computers to run simulations. Efforts are 
currently under way to convert these models to a more user friendly format, so that they 
may be used on smaller computers. The empirical models, however, are readily available, 
and take little computational time per simulation. One drawback is that the accuracy of 
the models is quite a bit less than that of the general circulation models. 

The history of the empirical earth atmospheric models dates back to the efforts of L. 

G. Jacchia in the late fifties. Once the early satellites such as Sputnik and Pioneer had 
been launched, immediate drag analysis was performed, and atmospheric models 
developed based on this analysis. The early models were crude, and only represented the 
regions where the satellites were orbiting. Little was understood of the variability of the 
environment and the density fluctuations being encountered. As more satellites, and a 
greater understanding of the sun's interaction with the earth's atmosphere was obtained, 
the accuracy of the atmospheric models increased. Figure 12 illustrates the developmental 
history of the various earth atmospheric drag models. (Marcos, 1993, p. 20) The 
mathematical development of the Lockheed-Densel or Jacchia 60 model, the Jacchia 71, 


and the MSIS 86 earth atmospheric models are described below. 
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1. Lockheed Densel 
The Lockheed Densel, or Jacchia 60 earth atmosphenc model was the first 

model to be implemented into the prediction scheme. As can be seen from Figure 12, it 
was one of the earliest atmospheric models contrived, and hence its accuracy is in great 
question. The Lockheed Densel model is actually a combination of two atmospheric 
models; the ARDC 1959 density model for low altitudes (h < 76 nautical miles), and the 
Jacchia 60 density formulation for h = 76 nautical miles. For the density below 76 
nautical miles, the density is obtained from the ARDC 1959 model which contains a table 
of density values as a function of altitude. The discussion on obtaining the density below 
76 nm will not be discussed in this paper. When the density at an altitude above 76 nm is 
desired, the Jacchia 60 formulation is used. The first requirement is to define the unit 
vector to the diurnal bulge using the solar position and the bulge lag angle. In order to 


accomplish this the solar longitude must be determined by the following equation, 


p=(2md/. .)-1.414+0.0335sin(274/_ ..} (15) 


where d represents the number of days since December 31, 1957. The unit vector to the 
sun is obtained from the following series of equations: 
&=cosA m., = sin A. COSE m= sin Asin € (16) 
where € is the obliquity of the ecliptic. (Lockheed, 1992, p. B-100) 
The unit vector to the diurnal bulge is then calculated in true of date coordinates 


by the following matrix 
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f,cos@—m, sin @ 
Up = C|m.cosé+fsin @ (17) 


Ns 

where C = J2000.0 to true to date transformation matrix, and 9 = 0.55 radians which 
equals the bulge lag angle. Two options are available to the user when using the solar flux | 
value of F10.7. The user may either specify a value, or one Is calculated using the 


following formula, 





07 = 1.5+.8cos274/ ., } (18) 


where once gain d is the number of days since December 31, 1957. This equation allows 
an approximation of the F10.7 value on any given day over the 1 1-year solar cycle period. — 
(Lockheed, 1992, p. B-103) The Jacchia 60 model divides the atmosphere into a series of 
three bands for density calculation. The first band is from 76 - 108 nautical miles. The 


equation used to calculate the density in this altitude band is given by 


fie on) es =)4 0.85F io {° =" 


«| 10+( 22 0+ cose) | (19) 
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p=71.18 
(p).. = density from ARDC 1959 model at 76 nm 


h = altitude in nm 


Fiy7 = solar flux measurement 


ReU; 





cosy = R = SV true of date position vector 


U; =diurnal bulge vector (equation 17). 


I) 


The second altitude band ranges from 108 - 378 nautical miles, and the density 
calculation is given by equation 20; 


p(t) = p(h)0.85F 107] 1.0+0.02375(e°"* -1.9)(1.0+cosg)"] (20) 


where poh) = k exp|(-b ~ab+ 6. 363e 0" | log.l 0] 


a=0.00368, 5 =15.738, and & is the conversion factor from slugs to kg/km: The third 
and final altitude band that is calculated in the Jacchia 60 model lies between 378 and 


1000 nautical miles. Equation 21 is used to calculate the density in this region. 
p(n) = (0.00504 Fie Yo. 125(1.0+ cos p)*(h? — 6.0x 10°)+6.0x 10° Jk (21) 
In the Jacchia 60 model it is assumed that the density is zero when the altitude is above 


1000 nautical miles. 


a. Mathematical Basis 
In order to provide a mathematical background for the above density equations, 
the partial derivative equations are illustrated below. The partial derivative of the 


calculated density with respect to position is described by equation 22; 


dh) anh) dh daplh) ag 
| ea sla ae (22) 


———SS re eee 


OR dh OR dp OR 
where kj is the conversion factor between nautical miles and kilometers. The partial 


derivative of the altitude with respect to R is best estimated by the following equation, 


Revl-e 


a ———— (23) 
l-—e° cos’ ¢’ 


R =position vector magnitude R. =earth equatorial radius 
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e = earth eccentricity o’ =geocentric latitude 


The differentiation yields equation 24; 


Oh RT 2 Vl—e*e* cosd’ | dcosd’ 
): OR 


OR R 3, 


(1 —e* cos’ ¢’ 


where 


: VX'+¥" dcose’ 


cos@’ = —————-,,_ and See 
? R OR R* cos@’ 





| XZ*,¥Z?,-Z(X?+¥*)| and 
X, Y, and Z represent the geocentric coordinates of the orbiting object. The partial of 
with respect to R is described by the following relationship, 


— — — 47 
we, — Us 
R3 R-=| (25) 


ag _ | 
OR sing 








The partial derivative of p(A) with respect to h and the partial of Ah) with 


respect to @ are different for each of the altitude bands. In the altitude band between 76 


and 108 nautical miles, the partial derivative of p(/) with respect to h and @ are given by 


the following equations: 


} a ate 
OAM) _ _P yy 2A I-L133F nf * =) 
oh 32 32 





h 


AB 3 
+——/]1+cos 26 
1274 | | a 


where Ah) is the density found from equation 19, p = 7.18, and the variables A, B, and C 


follow (Lockheed, 192, p. B-105): 
P 4/3 
A= (o)| | B= (Et) 0.857 (7 ) 
h 32 32 
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h-76 3 
c=} 1.04 pey \i-+eos0) | 


agli) 


== -345( 7S ls cos) sing (27) 
y 
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For the altitude band between 108 and 378 nautical miles, these two equations 


take on this form; 


doth) 
he 


= -(h)[a +.0.0305424 exp(-0.0048/) |é10 


+0.0002059125p.(h)Fi0.7{exp(0.0102h)|(1+cosg)’ (28) 


dpth) 
09 


= -0.0605625p.(h) Fio.7[exp(0.0102h)-1.9|(1+cosg) sing (29) 
In the highest altitude band between 378 and 1000 nautical miles, the equation 


take on the form (Lockheed, 1992, p. B-106) 


Axh)  — 8p(h) ,0.00189F 10.7 


l k 
= : ss [1+ cos]. (30) 


a = a +cos@) sin gh? -6.0x 10°) |A (31) 
These equations can be used to implement an atmosphere in an orbital 
prediction program, and in fact are currently used for that purpose. The Jacchia 60 
density model was one of the first empirical models on the market. With the launching of 
numerous satellites throughout the recent years, Jacchia et. al. have been able to expand 


on their empirical model. As can be seen from figure 12, several Jacchia models have been 


developed, each model building on its predecessor. 
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2. Jacchia-Roberts 71 

The next model to be discussed is the Jacchia-Roberts 71 atmospheric model. 
The Jacchia-Roberts 71 atmospheric model takes into account both the solar and the 
geomagnetic activity during the time in question. L. G. Jacchia defined two empirical 
profiles to represent temperature as a function of altitude and exospheric temperature. 
One profile was defined between the region of 90 to 125 kilometers, and the other above 
125 kilometers. Jacchia then used these temperature functions in the appropriate 
thermodynamic differential equations to obtain density as a function of altitude and 
exospheric temperature. (Lockheed, 1992, p. B-81) The Jacchia model as it stood was 
very cumbersome and required a great amount of data storage to hold the data required, 


so C. E. Roberts provided a method for evaluating the Jacchia model analytically. This 





lead to a faster and more manageable model. The following are the equations used in 
determining the atmospheric density as derived in the Jacchia-Roberts 71 atmospheric 
model. 

The first step in the model is to calculate the nighttime minimum global 
exospheric temperature for zero geomagnetic activity, 

Te = 379° +3°.24 Fo 1+1°.3| Fi0.7— Fio.7] (32) 

where fio 7= 81 day running average of the F10.7 centered on the day in question. Fo 7 
= solar flux measurement as obtained from the solar observatory at Ottawa, Canada. 
(Jacchia, 1970, p. 16) 

The value of the nighttime minimum exospheric temperature is then used tn 


calculating the uncorrected exospheric temperature as follows, 
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n= [4140.9 sin §+(cos** n-sin?* @)cos*® =| (33) 


where 
n=-|0-8 9=-l9+8 t= H -37°.0+6°.Osin( H +43°.0) , 
me TST 


6; = the sun's declination 


@ =the geodetic latitude of the satellite in true of date coordinates ( includes 


earth flattening) 


H = 180°. (SX, -5,%,) cos” = sa 7 (34) 
mSX2-S Xi | (5? 452)7( x7 + x2)" 


Fhe X variables are the components of the unit position vector of the satellite in true of 
date (TOD) coordinates, and the S variables are the components of the unit vector to the 
sun in TOD coordinates. (Lockheed, 1992, p. B-83) 

It has been found through analyzing the effect of geomagnetic activity on the 
atmosphere, that there is a lag of approximately 6 to 7 hours from a detection of a density 
change from the actual geomagnetic disturbance. In order to account for this lag, the 
value of Kp is obtained for a period of 6.7 hours prior to the integration time in question. 
It must be noted that the Kp value only exists at a three hour resolution. At this point the 
correction for the exospheric temperature is calculated using the following formulas, 

AT = 28°.0Kp +0°.03e*? (Z2 200 kilometers) 
AT==14°.0Kp+0°.02e*? = (Z< 200 kilometers) (35) 


The corrected exospheric temperature is then 


a] 


T~=1:+Alq~ 

and the inflection point temperature is given by the following formula (Jacchia, 1970, p. 
21) 

Tx = 371°.6673 + 0.518806 T~ — 294°. 3505e7 710" T= (36) 
These values for the temperatures and the satellite altitude are used in the calculation 
presented by Roberts for the atmospheric density. However, a number of corrections 
must be applied due to several atmospheric processes presented by Jacchia. These 
corrections will be presented before continuing with the Roberts calculations. 

One of the corrections deals with the direct effect of geomagnetic activity on the | 
density below 200 kilometers. This correction is calculated by the following relation, 

(Alogp)c =0.012Kp +1.2x10%e” (37) 

The next correction deals with the semi-annual density variation. This correction’ 

is calculated from the following relations, where Z 1s the altitude in kilometers. 
(Alog p)s1 = f(Z)g(t) (38) 

where 

F(Z) =(5.876 x 1077 Z?! +0.06328) e852 (39) 


g(t) = 0.02835+/0.3817+0.17829sin(2 wts1 + 4.137) ]sin(47ts1+4.259) (40) 


1.65 
ts = b+ o.ossa{| T+ Lsn(anes 6.035) (41) 
OW 
ant 71a 42 
O 365.2422 42) 


where JDysss is the number of Julian days since January 1, 1958. (Lockheed, 1992, p. B- 
84) 
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The next correction deals with the seasonal latitudinal vanations in the lower 
thermosphere. Equation 43 represents the general density variation, whereas equation 44 
represents the correction for Helium specifically. 


-0.0013(2-90)"] 


(Alogp)zr =0.014(Z- 90)e: sin(2 7+1.72)sin@|sin@| (43) 


Os|| . La dds _ 
a E = 0.25355 (44) 





(Alog p) He = 0.65% 











where € is the obliquity of the ecliptic. (Lockheed, 1992, p. B-84) 


Below 125 kilometers, Roberts uses the same temperature profile as Jacchia, 
dix z 
T(Z) =Tr+— > CuZ (45) 
35 n=0 


where 

di=Tx-To, To=183°.0K and the coefficients follow 

Co = ~89284375.0 C2 =—52687.5km~ Ca = -0.8hm7 

Ci = 3542400. 0km™' C3 = 340.5km~ (46) 
where 7; is the inflection point temperature at 125 kilometers calculated by equation 36. 
Roberts then substituted the temperature profile obtained from equation 45 into the 
barometric differential equation and integrated by partial fractions to obtain the following 
expression. Equation 47 represents the density found in the altitude band between 90 and 


100 kilometers. (Lockheed, 1992, p. B-85) 


p:(Z) 





= a ot exp( KF, } (47) 


where the "o" indicates the conditions at 90 kilometers. The mean molecular weight is 


calculated by the following equation, 


Zo 


M(Z) = Sn” (48) 


n=0 
where 
Ao = —435093.363387 A3 = 11.043387845hon~ 
At = 28275.564639 1 Am" As = -0.08958790995km— 
Az = —-765.33466108km~ As = 0.00038737586km~ 
Ae = —0.000000697444 kin 


These constants give a value of the mean molecular weight at 90 kilometers of 28.82678 
which is close to the sea level mean molecular weight of 28.960. (Lockheed, 1992, p. B- 
86) The density of the lower limit is assumed to be constant at p. = 3.46107 gm/cm’ | 


The constant & in equation 47 ts evaluated by the following formula 


sin 
RdiCa 


k= (49) 


where g is assumed to be 9.80665m/sec2 , the acceleration due to gravity at sea level. 


The functions /, and /, in equation 47 are determined from equations 50 and 51. 
P| Are hat Nel ZA Na a ee (50) 
' (90+R,} (90-7) \90-r,} (8100-180 + X*+Y° 


= Ps Ps -| Y(Z -90) 
F=(2-90) db] Bt en | ON 


The variables r, and r, are the two real roots and X and Y are the real and 





imaginary parts (Y>0), respectively, of the complex conjugate roots of the following 


quadratic, 
4 
P(Z)= >) Cn2" (52) 
n=0 


with the following coefficients 
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_ eee Go and Gas 


>= Isns4 
Crees C, 





for values of C, used in equation 45. The parameters p, in equations 50 and 51 are 


evaluated by the following relations (Lockheed, 1992, p. B-86) 


S(r) Lela) es) 
fae (r,) Ps; = Vy 











Sy 


Ps ={B,—nr,R?[B, + B(2X +r +r, —R,)|+W(r)p, }/X° 
-{nr,B,R, (X?+Y?)+W(r,)p, +n7,(R? -X?- Y?)p, jx" 
Pp = B, +B (2X +r, +r, -R,)—ps-2AX+R,)p,-(n+R,)p; -(7, + R,)p, 


B= 2 py") — Pp, 





where 
x's ~2rr,R,(R: +2XR,+X° ie) 
V =(R, +r, \(R,+r,)(R2 +2XR, + X72 +7?) 
Uae (7 Fux i (r/ DON, + X° +¥?\(r, -r,) 
2 ae Z 
W(r,) =nr,R,(R, eR ee] 
r 
and 
0s 
B.=a,+B, (n = 0, 1,.... 5) 
ly ~ 49 
5 
= 5/2: 
n=0 
with the coefficients 
Qt, = 31449025 16.672729 B, = -52864482.17910969 
Q, = -123774885.4832917 B, = -16632.50847336828 
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a, = 1816141.096520398 B. = -1.308252378125 


a, = —11403.31079489267 B, =0.0 
O, = 24.36498612105595 (Ee) 10, 
a, = 0.008957502869707995 p,=0.0 





Equation 47 is a valid equation below 100 kilometers where mixing 1s assumed 
to be predominant. However, above 100 kilometers, diffusive equilibrium is assumed, and 
it is necessary to substitute the profile given by equation 45 into the diffusion differential 
equation (one for each constituent of the atmosphere), and integrated by partial fractions. 
This was done by Roberts to yield the density for the altitude band between 100 and 125 


kilometers, given by equation 52. 


p,(Z) = >’ p,(Z) (52) 


It is computationally expensive to calculate the density at 100 kilometers at 
different exospheric temperatures by using equation 47. Thus values for the density at 100 
km were pre computed at several values of the exospheric temperature, and these values 
could then be extracted instead of using valuable computer time and memory. Next the 


atmospheric constituent mass densities are calculated by the following, , 


T(100) 
TZ) 





1+a, 
M 
p(2)=H(00) | | Fy exp(M kF,) (53) : 
S 


The index / corresponds to the values |] through 6 of the various atmospheric 


constituents of No, Ar, He, 02, O, and H, respectively. The constants found in equation 


: 
53 are the characteristics of these species and are tabulated in Table 3 of Lockheed 1992. : 
Hydrogen is not a significant constituent below 125 kilometers, so it is not included in the 


a2 


calculations. (Lockheed, 1992, p. B-89) The temperature at 100 kilometers is calculated 


from the following, 


4 
7(100) = 7, +Qd, where =35~5°C,(100)" =-0.94585589 (54) 


n=0 


and F, and F, are calculated by equations 55 and 56 
ie = Ze WiaZaAe) | 2 a2, : 27 2X7 + X* " (55) 
> (R.+100) (100-7, } (100-1, } (100? -200X +.X7+¥? 


I{Z-100) Me paqet] —_M(2=100)__ 
(Z+R,\(R, +100) y tan Fe ae (56) 








f,= 


where the parameters q, are calculated by 








U(r) 
7 —| 
BU) 
i 
Ee V 
q, = {l4n(R? - X?-¥*)g,+W(n,)q, +W(r4)q5 f° 


Bie 2 Nats) garnet Reagent Req; 

N= —244—- 93 — 92 
All of the variables in these equations have been defined earlier. 

For the region above 125 kilometers, it is still a valid assumption that diffusive 
equilibrium is dominant, but the temperature given by equation 45 is no longer valid. 
Jacchia defined the temperature profile of the upper altitude regions by the following 


empirical asymptotic function: 
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” ee | 
(7) eee ~1).a 0954 : - ZB asaoe-129" (57) 
(i wo ty ) 


In order to integrate the diffusion differential equations in closed form, Roberts replaced 


equation 57 with the following (Lockheed, 1992, p. B-90) 


(2) =1.-(T.- re — | F: "| (58) 


The parameter ¢ will be discussed later. 











Integration of equation 58 yields the following equation, which is valid for all of : 


the constituents except hydrogen. 








eh 1+@,+Y, 7 ap " 
{(Z)=p,(125)| = - 59 
p(2)=p,i2s( 4) "(E=2) 59) 
where 
2 = 
ye is it ly ae (60) 
RT, \ Ty — 1h 6481. 766 


At this point several corrections are made for the particular constituent densities 
due to seasonal changes. The value of the helium density that 1s calculated by using 
equation 59, must be corrected due to the seasonal latitudinal variation as given by 
equation 44. The specific form is presented below 

[Ps(Z)]. eceg = P3(Z)10°"8" (61) 
where / corresponds to the index of helium presented above. Above 500 kilometers, the 
concentration becomes significant, therefore it must be accounted for by the following 
146 +76 's 
p,(Z)= p,(s00} | (62) 


where the hydrogen density at 500 kilometers can be calculated from 
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p,(500) = = J 1731371 39-4- 5.508 Tse oR Ts00| (63) 


The temperature at 500 kilometers is calculated by using equation 58. The 
constituents are summed and the standard density above 125 kilometers is given by the 


following (Lockheed, 1992, p. B-94): 


p,(Z) =Ylz) (64) 


So far, the standard atmospheric density at anv given altitude has been calculated, 
but the densities calculated using equations 47, 52, or 64 must now be corrected for 
geomagnetic activity, the semi-annual variation and the seasonal latitudinal vanation by 
equations 37, 38, and 43 respectively. The effects of these variations can be summed 
logarithmically to obtain the following relation 

(Alogp),,, =(Alogp), +(Alogp).,+(Alogp),, (65) 
The final corrected density at any altitude can then be determined by 

p{Z) = p,10°8* (66) 

Due to the introduction of equation 58 in the place of Jacchia's equation (57) for 
temperature, the results of the Roberts portion of the model did not exactly concur with 
those that were observed by Jacchia. This effect was only in the upper portion of the 
model, whereas the lower altitude bands agreed exactly(as should be the case). This is 
where the @ parameter comes into play in equation 58. Values of the 2 parameter were 
calculated in order to obtain the best least squares fit of density versus altitude (125 - 
2500km) to the Jacchia tabulated data. The maximum deviation from the Jacchia density 
values is less that 6.7%. (Lockheed, 1992, p. B-94) The derivation and values of the 


parameter @ can be found in the Lockheed reference, and will not be discussed in this 


paper. 
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a. Mathematical Basis 
The following section contains the partial derivatives of the Jacchia- 
Roberts model. Starting from equation 66, the partial derivative of density with respect to 
local position can be found. This can then be used in orbital prediction schemes in order 
to better model the satellite's orbital trajectory. The following equation is a simple 
restatement of equation 66. 
AZ) = Ps(Z) AP con (67) 


The desired partial derivative then becomes 


dp HAp,) dp; 
SS = — +A —— 
ar °° OR Pc OR 





(68) 


The variation of the correction factor with respect to the satellite position is derived from 


equations 65 and 37 through 43. 


aa -_ j . g — + in(2 72 ~0.0013(Z-90) 








4 1-0.0026{(Z-90)}* sin gJsin ¢ Lae 90)|sin ¢| coset | (69) | 
OR OR | 


where 
f'(Z) = 0.002868 f(Z)+2.331(5.876 x 1077 )Z!3!e0 m28882 (70) 
The variation of altitude with respect to position, dZ/dR, is the same as that 
calculated in equation 24. The variation of geodetic latitude with respect to position is 


derived by differentiating the following equation 


| ne 


=f | 
ee O- [Oe 


(71) 


to obtain 
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st 

Xi EG 
ao _ sin2¢ » ee 
OR 2 X) + X3 


te! 


The variation of standard density with respect to position is calculated directly 


(72) 


from the barometric differential equation for altitudes below 100 kilometers 


Ops all pel Yin, Sj 1) (73) 





OR OR TOR 


and from the diffusion differential equation above 100 kilometers 


Ap, ] p’g,R, | aZ oT 
[ee | ate 74 
ae a has at (PstOsPsdaer (74) 


where 


6 
=) P.M, (75) 
1=] 


The partial derivatives of the temperature with respect to position are derived by 


differentiating equation 45 for altitudes less than 125 kilometers 


or T-T,{ oT, aT. - \oaZ 
=—— — : 12 76 
OR d, Zi a (Sey \g 7) 


or equation 58 for altitudes above 125 kilometers 


oe ele | pee -|72\ =) 
OR OR \T.-T, of,} \R.+z\\35 
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ae - R,+Zy, |{ € \oz 
teehee 


and finally, the derivatives of 7, and 7, with respect to position are derived by 


differentiating equations 36 and 33 respectively, 





= = 0.0518806 +(294.3505)(0.0021622)e "627 (78) 


os 


or, = (0 RAS ie aos —cos*” =\e 
a 2H 'OR 





-2.2cos'* nsin Ncos”® ESA 3 cos" n- sin?? 8) cos” sins (79) 
where 

(ait = eee 

OR 2|¢-6,| OR 

00 _ 1 o+6, db 

OR 2|¢+65,| OR 
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3. MSIS 86 (Mass Spectrometer and Incoherent Scatter 86) 
The Lockheed Densel model and the Jacchia 71 model, both developed by 

Luigi Jacchia, use drag analysis data from earth orbiting satellites to derive the densities of 
the various altitude bands associated within the models. A. E. Hedin chose a different 
approach to the problem of modeling the atmosphere. Hedin's approach was to use actual 
density data retrieved from orbiting instruments aboard several satellites and combine this 
data with measurements taken by ground based incoherent scatter radar stations. Hedin, 
et. al., began the original MSIS model using the total densities inferred from Jacchia's 
models, however, it was found that the measurements of the temperature found by the 
incoherent scatter method differed significantly. (Hedin, 1972, p. 2139) These differences 
were noted in the time of daily maximum, and the amplitude of the daily and annual 
variations of the density. Hedin used the data that was obtained from the mass 
spectrometers on board several orbiting satellites to confirm this. The measurements from 
the satellites as well as the ground based incoherent scatter measurements provide 
information at different altitudes, latitudes, longitudes, solar activities , and seasons. 

The MSIS 86 model is the culmination of years of research and development as 
can be seen in figure 1. The formulation of the MSIS models is based on a spherical 


harmonic expansion of exospheric temperature and effective lower boundary densities 
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which uses local solar time and geographic latitude as the independent variables. This 
expansion is a special case of a more general expansion based on longitude and latitude, 
where the coefficients of the expansion are represented by a Fourier series in universal 
time. (Hedin, 1979, p. 2) 
This leads to the following expansion 
AG = x y yP alae cos(ma’t +nA) +b" sin(mo’t +n)| (80) 
m=0 1=0 n=-/ 
where 

P,, = Legendres associated functions in geographic latitude 

f = universal time in seconds 

cw’ = 27/86400s"' 

dX. = geographic longitude in radians 

This particular expansion was used for the initial MSIS models. With the 
launching of several more satellites equipped with mass spectrometers, and the addition of 
several more incoherent scatter radars, the model evolved into the MSIS83, and eventually 
the MSIS 86. 

The main emphasis now is the formulation of the MSIS 86 model. Hedin 
explains that the models uses a Bates temperature profile as a function of geopotential 
height for the upper thermosphere, and an inverse polynomial in geopotential height for 
the lower thermosphere. (Hedin, 1987, p. 4649) The exospheric temperature and other 
key quantities are expressed as functions of geographical and solar/magnetic parameters. 
The temperature profiles allow for the exact integration of the hydrostatic equation for a 
constant mass to determine the density profile based on a density specified at 120 
kilometers as a function of geographical latitude and solar/magnetic parameters. The 


MSIS 83 model used the expansion formula (Equation 80) to model the variations due to 
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local time, latitude, longitude, universal time, F10.7, and Ap. The MSIS 86 model 
enhances the MSIS 83 model by adding terms to express hemispherical and seasonal 
differences in the polar regions and local time variations in the magnetic activity effect. 
The following equations are the complete formulation of the MSIS 86 
atmospheric model. The first step in the formulation is to obtain an expression for the 
temperature profile. 
AT —(T. ~T) exp! o(2.2))| Where Z > iz (81) 
1(Z)= Mu T, x? +Tox*+Tx°) whereZ<Z, (82) 


the matching temperature and temperature gradient at Z, equals 


, =1(Z,) =, -(L, -T exp 4") (83) 

I, =1(Z,)=(L.-T)o(R,+2,)/(R,+2,)) 84) 

i —- 2) 0/7? -3.1III/T, -Y7)+7.1111 11/7, -/7,) (85) 

Z,,2.)T\27,)-(/T, - Wh) - 2% (86) 

a. -T, (87) 

r= —(Z,,Z)-&Z,.Z,)|/6(Z,.Z,) (88) 

aan, = 7) (89) 

(Z,Z,)=(Z-Z,)(R, +Z,)/(R, +Z) (90) 

{(Z,Z,)=(Z-Z,R, +Z,)/(R, +Z) (91) 
=1, [(T,-7) f= T[1+G(1)] 
T, =7, [1+G(L)] I, = T[1+G(Z)] 
Tl, = T.[14+G(Z)] Z, =Z,[1+G(Z)] 


where (all temperatures in Kelvin) 


T =ambient temperature [, = average mesopause temperature 


4] 


[, =temperature at Z, 7, = average temperature gradient at Z, 


, 


[, =temperature at (Z, + Z,)/2 Z = altitude in km 
T. = average exospheric temperature Z, = (2onm 
7, = average temperature at Z, 


Z,, =altitude of temperature profile junction; 117.2 km 
Z, = average mesopause height R, = 6356.77 km 


Hedin now explains that the density profile is a combination of diffusive and 
mixing profiles multiplied by one or more factors C,...C, to account for chemistry or the 


dynamics flow effects. (Hedin, 1987, p. 4638) 


VA 


n(Z,M)= [n,(Z, M)* +n_(Z, M)"| C,(Z)... ca)" Kem 


A=M,/(M,-M) (93) 
where 
n= ambient density M,, = 28 
n, = diffusive profile M, = 28.95 
= mixing profile M = molecular weight of gas species 


The Following equations represent the diffusive profile and the mixing profile. 
The diffusive profile is given by the following 
n,(Z,M) =n,D(Z, M)[1(Z,)/7(Z)}" (94) 
DZ) = sD a7 ee (95) 
D(Z,M) =D, (Z,,M)exp err ipereleyserel NP (7 7) (96) 


iD AZ M =e (TA) ]" ex pes Z Zu) (97) 
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T. =temperature gradient at Z, in K/km ip. = average mesopause shape factor 





a= Mg,/(oR,T.) Sa = g,/(1+Z,/R,) 


y, = Mg,&(Z,,Z,)/R, g, =,/(1+Z,/R,) 


n=, exp (4) (98) 
where 
n, = the average density at Z, 
g, =9.80665x 10° km/s 
R, =8.314x10> g km2/mols s* deg 


and « is the thermal diffusion coefficient of -0.4 for He and zero for the rest of the 


constituents. 


The mixing profile is given by the following 


y(Z,M) =n, | D(Z,,M)/D(Z,,Mo)[7(Z,)/1(Z,)] (Z,Mo)[7(2,)/7(Z)] (99) 


Hedin (1987), goes into a procedure to simulate the chemistry and dynamic flow 


effects of the turbopause which provides a specified mixing ratio with respect to nitrogen. 
C(Z) = expLR, /[1 +exp(22i/hi } (100) 
R, = log|Qu,,(Z,,28)/n,,(Z,,.4)| (101) 
where (2 is the mixing ratio relative to N2, Z] 1s the altitude at which log Cy 1s R/2, 


and Hj is the scale height for this correction. There is a peak in the mixing ratio in the 


lower thermosphere due to O, N, and H. This peak 1s modeled by the following 
C,(Z) = exp| Ry /[1 + exp(Z- Z,)/H, ]} (102) 
where R> Is the density correction parameter, Z> 1s the altitude where log)9 density 


correction is R»/2, and H> 1s the scale height of the correction. As mentioned above, the 
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MSIS 86 atmospheric model attempts to model all the known causes of the variations of 
density. This is done by the expansion function (Equation 80). The following is the 


expansion function for the MSIS 86 model quantities: 


Time independent 

G=a5 a lo sta wiee 
Solar activity 

+F AF + F22( AF) + fOlAF + fo (AF) + f2'P,,AF 
Symmetrical annual 

+c! cosQ,(t, - 13) 


Symmetrical semiannual 





+(c3, +3, Po) cos22,(t, ~ 163) 


Asymmetrical annual (seasonal) 


+(cl,P, +cl,Py)F, cos2Q,(1, - 17!) 


Asymmetrical semiannual 


Cg cos2Q,(t, = (3) 
Diurnal 
+/4,,P, +ayP, +4yP,, +04, P, cosQ,(t, 13) |F, coset 
+/,P, PO a 1 One ten cos, (1, mG IF, sin OT 
Semidiurnal | 
+a. Py tay Fy, +(c3, ree P,,)cosQ, (1, =the IF, cos2@T 


+[b,, Poy +4 Pay +(dj, Py +5, P,, )oosQ, (1, - 13) |F, sin 20 
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Terdiurnal 


+|a55P, +(ci,Ps +¢4,P,,)cosQ,(t, — to) [Fa COsS3@T 
+52 +(di,P, + do, P,,)cosQ,(t, — if F, sin 3@T 


Magnetic activity 


| fey + he Pay + Hig Pag + (esa Pro + 5p Poo + ksp Pop) COSQ, (1, — Hi) 
+k, +45,By +k, B,) cosa t- 1, )|A4 


Longitudinal 


0 0 0 0 0 0 ] l | 
+a, P, +a),P, +45,P, +a,P, +43, P, +a, P, +(a A, +45 P, )cosQ,(t, -13)| 


x(1+ F°AF) cosa 


+/B9,P,, +9, Py +O5,P, +5 ,P, +b,P, +55,P,, +(d4P, +85 P, )cosQ, (1, -3)| 
x(1+ FEAF) sin A 
UT 
+(a!,P. +a!,Py+a',Po)(1+ Fg AF )(1+ 72!P,)(1+ 72! cosQ,(¢, -12)) 
xcosa’(t—12')+(al,P, +a3,P.,) cos o'(t —121)4 22] 
UT/longitude/magnetic activity 
HA Py +k Py + ke P (+14? Ro JAA cos A ~ At?) 
+k! P cosQ,(t, — 1 )AA cos A— A] 


+{ kein Py +k Py +h) Pay |AA cos’(1 — th!) (103) 
Fo=1+FAF + fol aF + fe (AF) (104) 
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Fy = 14+ FiAF + fa aF + fo (AF) (105) 


where 


AFy="F eae AF = F.,, -150 F10.7 is 
the solar flux measurement for the previous day, F.,, is the 81day average of the solar flux 


measurement centered on the day in question. 


Pare the non normalized Legendre-associated functions equal to the following 
(1 = /2'n\ aren /ate™*™)(x -1)' 
with x equal to the cosine of the geographic latitude, 
Q., =27/365 w@=27/24 w’ = 27/86400, Tt is local time in hours, t is 
universal time in seconds, tq is the day count in year, and A is the geographic longitude 


with eastward being positive. Hedin explains that there are two choices for magnetic 


activity, 
AA =(Ap -4)+(ki, - 1)\ Ap -4+|exp(—z¢.( Ap =a 1/ ke} (106) 


where Ap is the daily magnetic index, or there is another method illustrated that 
establishes the magnetic activity that uses the average values over an extended period of 


time. (Hedin, 1987, p. 4661) 
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II. ATMOSPHERIC MODEL EVALUATION 


The evaluation of atmospheric drag models is a relatively difficult procedure owing 
to the fact that there is no real comparison tool to be used. Until enough satellites can be 
placed into various orbits and each of these satellites carries density measuring equipment, 
the density of a particular orbit cannot truly be obtained. As mentioned previously, the 
modeling of the earth's atmosphere is the most inaccurate of all the perturbations 
encountered in orbital motion analysis. The intent of this project was to obtain known 
vectors of a particular satellite program, and run a comparison of an orbital prediction 
scheme using each of three atmospheric models for two different altitudes. In order to 
make the comparison valid, consistency of the prediction scheme had to be verified. 

The vectors being used in the analysis were obtained from the Air Force Satellite 
Tracking Center on two satellites, one at 650 kilometers and one at approximately 800 
kilometers, both of which have a high inclination. The vectors, position and velocity 
represented in True of Date Cartesian coordinates, are obtained by triangulating the 
satellite position and predicting out to 20:00:00.00 local time. It would have been better 
for the evaluation, if the actual observed position of the satellites could have been 
obtained, but due to time constraints and logistical problems this could not be achieved. 
The accuracy of the position and velocity vectors is advertised by the STC to be within 
2000 ft in the X, Y, and Z directions. It was decided at the time, that this accuracy was 
acceptable for the evaluation. 

In order to provide consistency, the prediction scheme developed had to be the same 
as that used by the STC. This proved to be quite a troublesome problem. The STC uses a 


prediction scheme which models most if not all of the known perturbations. Depending on 
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what is requested, certain perturbations can be temporarily removed from the program. 
The vectors obtained from the STC contained only those perturbations which were 
inherent to the prediction scheme developed, with one exception. As mentioned 
previously, the STC uses the WGS-84 geopotential model. A great deal of effort was 
placed on placing the same geopotential model in the developed prediction scheme, but 
several problems soon became apparent. The conversion of True of Date coordinates into 
the WGS-84 frame of reference is not a trivial process. An attempt was made to 
incorporate the full 41 by 41 geopotential matrix in the scheme, but error build up and 
computational time increased greatly. It was decided to incorporate a truncated version of 
the WGS-84 geopotential model. A six by six matrix representing the first six zonals of 
the model was eventually used to prevent error build up. Future work on this project 
would have to include the conversion of the TOD coordinates to the WGS frame of 
reference and incorporation of the full 41 by 41 matrix. 

Another way that consistency was achieved was to verify and use astrodynamical 
constants consistent with AIAA standards. Depending on which reference one uses to 
obtain values such as the earth's gravitational parameter, earth's equatorial radius, earth 
flattening constant and others, these values would all vary by some small amount. It was 
found that depending on which constant was used, the error would be different. It was 
decided to use the AIAA published values whenever possible. These values were placed 
in a subroutine of the prediction scheme to be used as necessary. 

Once it had been determined that the developed prediction scheme was an equal 
representation of the STC's prediction scheme, or as close as it could be gotten in the time 
allowed, the model evaluation was begun. The vectors obtained from the STC spanned 
the time frame from 21 August 1993 to 22 September 1993 in 24 hour forecasts. That is, 


the position and velocity were reported for 20:00:00.00 of each day. These vectors were 
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then programmed into 60 satellite files, 30 for the 650 kilometer satellite and 30 for the 
800 kilometer satellite. Each of the satellite files contains not only the position and 
velocity vectors, but also the environmental conditions of F10.7 and Ap for that day, plus 
the 81 day average F10.7 as reported by the National Geophysical Data Center. The 
vectors were propagated over a 24 hour period and compared to the next day's vector. 
Since the baseline vectors themselves contained error, it was decided to compare the 
results by using the specific mechanical energy, a conserved quantity, of each of the 
satellites over the 30 day period. By converting the position and velocity of the STC 
vectors and those obtained from the developed prediction scheme into specific mechanical 


energy using equation 107, a comparison could be made. 
E= mn ba (107) 


A relative comparison of the position, velocity, right ascension, declination, and radial 
distance was then conducted in order to determine which model provided the most 
accurate results, the object being to meet or beat the 2000 foot accuracy advertised by the 
STC prediction scheme. 

As mentioned, the evaluation procedure consisted of propagating the baseline vector 
for a 24 hour integration period for each of the satellites, then switching atmospheres and 
running the routine again. It must be noted that the Jacchia 60 and Jacchia 71 models 
were run with relatively little difficulty. The MSIS 86 atmospheric model proved to be 
quite the opposite. Due to an unexpected hardware failure late in the project, and the fact 
that the atmospheric model itself is quite difficult to integrate into the prediction scheme, 
the evaluation of this model could not be completed in time. Work is currently under way 


to correct this situation. 
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A. JACCHIA 60 MODEL EVALUATION 

The first evaluation was conducted on the Jacchia 60 earth atmospheric model. The 
first procedure was to determine the values for the solar flux and geomagnetic activity 
encountered during the time period between August 21 to September 22, 1993. Figure 13 
shows the reported values of F10.7 and Ap for this time period. As compared to Figure 6, — 
it can be seen that the solar flux values are on the low side of the solar cycle, and the 
geomagnetic activity 1s relatively low. This allows for an easier comparison, because a 
perturbed atmosphere would have made some of the data points obtained during an active 
period out of range of those obtained during a quiet period. 

The first evaluation was conducted on the satellite at 650 kilometers. Figure 14 


illustrates the specific mechanical energy of the STC reported vectors versus the predicted 


vectors. It can be seen that the predicted vectors are extremely close to that of the STC 





reported values. The two error bars located on the graph illustrate the 2000 ft error range 
advertised by the STC prediction scheme. In one case it was found that the STC obtained 
vector on Julian date 242 was actually reported incorrectly. For the most part, however, 
the predicted values showed only small variations with the baseline vectors. In-track 
errors fell within the range of 1700 to 9000 feet in radial distance, and 0.01 to 0.15 
degrees in declination. Cross-track errors also fell within the range of 0.01 to 0.15 
degrees. Satellite velocities had some small error in the range of 0.2 to ten feet per 
second. 

The same procedure was then conducted on the 800 kilometer satellite vectors. 
Once again the prediction scheme provided vectors very similar to those obtained from the 
STC. Figure 15 illustrates the results. It must be noted at this point, that the accuracy of 
the prediction scheme suffered at the higher altitude. This is due to the atmospheric 


model. As mentioned previously, accuracy of the atmospheric model deteriorates the 
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higher the altitude. In-track errors were found to range from 2500 to 9000 feet in radial 
distance, and 0.01 to 0.2 degrees in declination. Cross-track errors were in the range of 


0.01 to 0.8 degrees. Velocity was in the range of one to eight feet per second. 


B. JACCHIA 71 MODEL EVALUATION 

Once again the procedure was run, this time with the Jacchia 71 model inserted into the 
prediction scheme. Figure 16 illustrates the results of the prediction runs on the 650 kilometer 
satellite. The prediction scheme provided very little variation with the baseline vectors. The in- 
track errors were 250 to 8000 feet in radial distance and 0.001 to 0.14 degrees in declination. 
Cross-track errors were also 0.001 to 0.14 degrees in mght ascension. Errors in velocity predicted 
were 0.1 to ten feet per second. The results of this atmosphere at 650 kilometers were much more 
accurate than those of the Jacchia 60 model. In order to confirm this, the prediction scheme was 
run again with the vectors of the 850 kilometer satellite. In-track errors were in the range of 200 to 
8000 feet in radial range and 0.001 to 0.10 degrees in declination. Cross-track errors were also 
comparable with the error range being 0.001 to 0.10 degrees of right ascension. This confirmed 
the fact that the Jacchia 71 model provided better accuracy at both altitudes. This is due to the fact 
that the Jacchia 71 model provides a more accurate modeling of the polar regions, whereas the 
Jacchia 60 does not model the polar regions as well. The modeling of the polar regions is 


extremely important for the satellite program being studied due to its high inclination. 
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IV. SUMMARY 


This thesis has attempted to provide the underlying principles used in orbital 
prediction and a comparison of atmospheric models. The prediction scheme that has been 
created was used to compare the Jacchia 60 and Jacchia 71 earth atmospheric models. It 
was found that the Jacchia 71 model provided much more accurate results than the Jacchia 
60. This would stand to reason due to the fact that the J71 model had a bigger data base 
with which to work from. The J71 model provides density information that is dependent 


both on solar activity as well as geomagnetic activity, whereas the J60 only relies on solar 





‘ 


activity. The J71 also models the polar regions, which is essential to the satellite program 
being investigated. The MSIS 86 earth atmospheric model could not be compared due to 
time constraints, but is currently being integrated into the prediction scheme. 

Several improvements can be made to the prediction scheme in order to improve 
accuracy. As discussed earlier, the geopotential model needs to be incorporated. The 
True of Date coordinates need to be transferred into the WGS-84 frame of reference, and 
the complete 41 by 41 geopotential matrix incorporated. Doing this would increase 
computational time, but accuracy would greatly improve. Another improvement in the 
prediction scheme would be to vary the B-factor as the satellite cross-sectional area 
changes through its orbital pass. Currently the B-factor is held at a constant value. By 
changing the B-factor as a function of latitude or declination, error can be reduced. 

Follow on work on this project would include inserting other atmospheric models 
for comparison, such as the MSIS 86 and MSIS 90. Before doing this, however, it would 
be better if the satellite vectors being used were actual observations and not predictions 
themselves. This would provide a better realization of the accuracy of the program, and 


hence the atmospheric models. Another improvement in the program's accuracy would be 





to convert from True of Date Cartesian coordinates to normalized spherical coordinates 


and then integrate. Doing this would greatly decrease the error build up. As the 
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prediction scheme stands now, the error build-up causes accuracy problems if the vectors 
are integrated for more than a twenty four hour period. 

Orbital prediction is an essential tool in today's space age. In the case of the satellite 
program of interest, if the accuracy of the orbital prediction can be improved, then the 
mission analysis can be improved. The interesting point of this project, was that the B- 
factor was no longer being used as the error catch-all. The B-factor for any given satellite 
can be correctly used, and the position and velocity as well as any other of the satellite 
ephemirides can be calculated. By incorporating more subroutines that model other 


perturbing forces, this prediction scheme can be used for other satellite programs. 
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APPENDIX A 
FIGURES 
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Figure 1. Atmosphere Breakdown 
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Figure 2. Temperature vs. Altitude 
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Figure 3. Temperature Profile of Earth's Atmosphere 
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Figure 4. Height Variation of Number Density 
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Figure 5. Diurnal Vanation in Density 
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Figure 7. Density vs. Altitude for Various F10.7 Values 
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Figure 8. Solar Radiation Spectrum 
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Figure 9. Geomagnetic Storm Recording 
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Figure 10. Density vs. Altitude for Geomagnetic Storm 
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Figure 11. Comparison of Perturbation Magnitudes 
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Figure 12. Atmospheric Model History 


65 


Measured Values 


Daily Measured Values of F10.7cm and Ap 


100 


60 OM ONGECHONONDE O10 oh Ono cor oO Oo Doo FDO GMO hb OOOH DODO OOS 6 
: 
. 
= Ap 
= 

* . . 
4 ; P ,SokSCSUenoOodoonmoOomoco noob oo eo namoea oo aco 

: A : . ' : 

i 5 : A ‘ : 

; 
20 ter eo eT re er OHS AR A Doe aD eee 





235 240 245 250 255 260 265 
Julian Date 


Figure 13. Observed values of F10.7 and Ap for 21 August to 22 September 1993 
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Figure 14. Jacchia 60 Evaluation at 650 km 
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Figure 15. Jacchia 60 Evaluation at 800 km 
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Figure 16. Jacchia 71 Evaluation at 650 km 
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Figure 17. Jacchia 71 Evaluation at 800 km 
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APPENDIX B 
TABLES 


TABLE 1: PERTURBATION ACCELERATIONS 
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TABLE 2: TYPICAL BALLISTIC COEFFICIENTS FOR LEO SATELLITES 
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